Boyé et al. 2019 Diversity and distributions

This is the code used to make the analyses of our article entitled Trait-based approach to monitoring marine benthic data along 500 km of coastline published in Diversity and distributions in 2019.

The article is open-access and can be found on the webpage of the journal at doi:10.1111/ddi.12987

Note that all figures that do not appear in this document were made on Inkscape. This include :

Some of the figures were also combined on Inkscape, including :

Please report any bug that you could find (better late than never!). Hopefully, there are none 🤞 🙏 🤞

Load home-made functions


Format the data

Community data

Data selected :

  • Spring
  • 2007 / 2010 / 2013
  • Infauna
# Select data and format the df
#------------------------------

faune_sel <- faune %>% 
  # Select Spring data, years 2007 / 2010 / 2013, and remove the epifauna ("Haveneau" = dip nets)
  filter(saison=="Printemps" & annee %in% c(2007, 2010, 2013) & method_echant != "Haveneau") %>%
  # Join species classification
  left_join(.,sp_classif,by=c("species"="species_init")) %>%
  # Modify the name of the habitats
  mutate(mini_theme=recode(theme, "Subtidal Meuble"="SM","Intertidal Meuble"="IM","Herbiers Intertidaux"="HI","Bancs de Maërl"="MA")) %>%
  # Change the name of the two maerl beds within the Bay of Brest: Rozegat & Keraliou
  mutate(site = case_when(
    str_detect(point_name,"Rozegat") ~ "Rade de Brest - Rozegat",
    str_detect(point_name,"Keraliou") ~ "Rade de Brest - Keraliou",
    TRUE ~ as.character(site)
  ))

# Aggregate data at the site level and format the site x species matrix
#----------------------------------------------------------------------

# Clean data
faune_cl <- faune_sel %>%
  # Aggregate the different entries that a single species may have in a sample
  group_by(theme,mini_theme,site,method_echant,annee,Class,species) %>%
  summarise(abondance=sum(abondance)) %>%
  ungroup() %>%
  # Remove the site of Keraliou with incomplete data
  filter(site!="Rade de Brest - Keraliou") %>%
  # Clean site names and order the sites and habitats
  mutate(site=recode(site,"Saint-Brieuc"="Baie de Saint-Brieuc")) %>%
  mutate(site=ifelse(site=="Baie de Saint-Brieuc" & theme=="Bancs de Maërl","Arcouest",site)) %>%
  mutate(theme=recode(theme, "Bancs de Maërl" = "Subtidal maerl beds", "Subtidal Meuble" = "Subtidal bare sediments", "Herbiers Intertidaux" = "Intertidal seagrass beds", "Intertidal Meuble" = "Intertidal bare sediments")) %>%
    mutate(theme = factor(theme,levels=c("Subtidal maerl beds","Subtidal bare sediments","Intertidal seagrass beds","Intertidal bare sediments"), ordered=TRUE)) %>%
  mutate(site=factor(site,levels=c("Baie du Mont Saint-Michel","Saint-Benoit","Saint-Malo","Saint-Briac","Saint-Cast","Baie de Saint-Brieuc","Arcouest","Sept-Iles","Lannion","Saint-Efflam","Pierre Noire","Morlaix","Callot","Sainte-Marguerite","Molene","Blancs Sablons","Iroise","Camaret","Roscanvel","Rade de Brest", "Rade de Brest - Rozegat","Rade de Brest - Larmor","Baie de Douarnenez Nord","Plage de l'Aber","Baie de Douarnenez Sud","Audierne","Glenan","Concarneau","Trevignon","Gavres","Erdeven","Lorient Etel","Plouharnel","Quiberon","Belle-ile","Meaban","Kerjouanno","Arradon","Vilaine Large Sud","Vilaine Large Nord","Vilaine Cote","Damgan"),ordered=T)) %>%
  # Transform character variables to factor type
  mutate_if(is.character, as.factor)

# Create the community matrix for the whole community
#----------------------------------------------------

faune_tot <- faune_cl %>%
  # Remove the Class variable
  select(-Class) %>%
  # Create site x species matrix by putting 0 for absent species
  spread(species,abondance,fill=0) %>%
  # Remove the precision for the azoique samples
  select(-contains("azoïque")) %>%
  as.data.frame()

# Identifiers
id_tot <- faune_tot %>%
  select(theme,mini_theme,site,method_echant,annee) 

# Community matrix
com_tot <- faune_tot %>%
  set_rownames(paste(.$mini_theme,.$site,.$annee,sep="_")) %>%
  select(-theme,-mini_theme,-site,-method_echant,-annee)

# Create the community matrix for the polychaete only
#-----------------------------------------------------

faune_pol <- faune_cl %>%
  # Select polychaetes
  filter(Class=="Polychaeta") %>%
  # Remove the Class variable
  select(-Class) %>%
  # Create site x species matrix by putting 0 for absent species
  spread(species,abondance,fill=0) %>%
  # Remove the precision for the azoique samples
  select(-contains("azoïque")) %>%
  ungroup()

# Identifiers
id_pol <- faune_pol %>%
  select(theme,mini_theme,site,method_echant,annee)

# Community matrix
com_pol <- faune_pol %>%
  as.data.frame %>%
  set_rownames(paste(.$mini_theme,.$site,.$annee,sep="_")) %>%
  select(-theme,-mini_theme,-site,-method_echant,-annee)

# Separate intertidal and subtidal data (different sampling tools)
#-----------------------------------------------------------------

# Intertidal polychaetes
com_pol_inter <- com_pol %>%
  rownames_to_column('id') %>%
  filter(id_pol$method_echant == "Carotte") %>%
  column_to_rownames('id')

# Remove empty species
com_pol_inter <- com_pol_inter[,colSums(com_pol_inter)>0]

# Subtidal polychaetes
com_pol_sub <- com_pol %>%
  rownames_to_column('id') %>%
  filter(id_pol$method_echant == "Benne") %>%
  column_to_rownames('id')

# Remove empty species
com_pol_sub <- com_pol_sub[,colSums(com_pol_sub)>0]

Trait data

i) Missing data exploration (Figure not included in the article)

  • Total number of missing data
Status Number of data
Available 9923
Missing 461
  • Modalities with missing data
Modality Number of missing data % of missing data
Short_life_span 115 48.728814
Medium_life_span 115 48.728814
Long_life_span 115 48.728814
Iteroparous 21 8.898305
Semelparous 21 8.898305
Dev_asex 17 7.203390
Dev_direct 17 7.203390
Dev_plankto 17 7.203390
Dev_lecitho 17 7.203390
Hermaphrodite 3 1.271186
Gonochoric 3 1.271186
  • 50% of missing data for Life span \(\Rightarrow\) Trait removed from the analysis
  • 9% of missing data for Reproduction frequency \(\Rightarrow\) Missing data were imputed
  • Distribution of missing data

ii) Format the fuzzy-coded trait matrix

# Format the fuzzy-coded trait matrix
#-------------------------------------
trait_fuz <- trait_mat %>% 
  # Remove the empty bioturbation modality and the life span trait with lots of missing data
  select(-Bioturb_R,-ends_with("life_span")) %>%
  # Select only the species in the community matrix
  filter(species %in% colnames(com_pol)) %>%
  # Order the species as in the community matrix
  slice(match(colnames(com_pol),species)) %>%
  # Define row names and format the final trait matrix
  as.data.frame() %>%
  set_rownames(.$species) %>%
  select(-Kingdom,-Phylum,-Class,-Order,-Family,-Genus,-Species,-species)

# Recode the matrix to give the same weights to all traits 
#----------------------------------------------------------

# Create the function to convert the fuzzy coded matrix (going from 0 to 4 for each modality) to a relative coding in which all modality of a given trait sum to 1 for each species
transfo_traits <- function(x) {
  if (is.na(sum(x))) { # If we have NA's, gives back NA's
    rep(NA, length(x))
  } else if (sum(x) == 0) { #  If we have only 0's, gives back only 0's
    numeric(length(x))
  } else { # Else, we get the modality's score and divide it by the sum of the scores of the species for all modalities of this trait
    x / sum(x)
  }
}

# Prepare the trait list to apply the function
trait_list <- trait_list %>%
  # Remove the empty bioturbation modality and the life span trait with lots of missing data
  filter(Accronyme != "Bioturb_R" & Trait != "Life span") %>%
  # Reorder the list according to the order of the modalities in trait_fuz
  mutate(Accronyme=str_replace(Accronyme,"-",".")) %>%
  slice(match(colnames(trait_fuz),Accronyme)) %>%
  # Reorder the trait factor levels
  mutate(Trait = factor(Trait, levels=unique(Trait))) %>%
  droplevels()

# Compute the number of traits
n_trait <- nlevels(trait_list$Trait)

# Define the start and end point of each trait
n_mod <- trait_list %>%
  # For that we compute the number of modalities for each trait
  group_by(Trait) %>%
  summarise(n=n()) %>%
  pull(n) %>%
  # And we do the cumulative sum
  cumsum(.) %>%
  c(0,.)

# Apply the function
trait_prop <- trait_fuz 
trait_prop[,] <- NA

for(i in 1:n_trait){
  # Select all the modalities of the trait i
  tmp <- trait_fuz[,seq(from=n_mod[i]+1, to=n_mod[i+1], by=1)]
  # Apply the function to this trait
  tmp <- t(apply(tmp, 1, transfo_traits))
  tmp <- as.data.frame(tmp)
  # Save
  trait_prop[,seq(from=n_mod[i]+1, to=n_mod[i+1], by=1)] <- tmp
}

summary(trait_prop) %>%
  kable("html") %>%
  kable_styling(full_width = F, position="left") %>%
  column_spec(column = n_mod+1, border_right = TRUE, include_thead = TRUE)
Size_inf2 Size_2.5 Size_5.10 Size_10.50 Size_50.100 Size_100.200 Size_sup200 SSDF SDF ASF PSF Grazer Pred Scav Parasitic Microphagous Macrophagous Infaunal Epibenthic Tube_dweller Burrower Crawler Swimmer Attached Mob_0 Mob_inf10 Mob_10.100 Mob_100.1000 Bioturb_N Bioturb_S Bioturb_B Bioturb_UC Bioturb_DC Hermaphrodite Gonochoric Dev_asex Dev_direct Dev_plankto Dev_lecitho Iteroparous Semelparous
Min. :0.000000 Min. :0.00000 Min. :0.00000 Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.000 Min. :0.00000 Min. :0.000 Min. :0.00000 Min. :0.00000 Min. :0.00000 Min. :0.0000 Min. :0.0000 Min. :0.000000 Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.00000 Min. :0.000000 Min. :0.000 Min. :0.0000 Min. :0.0000 Min. :0.00000 Min. :0.00000 Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.00000 Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.000 Min. :0.0000 Min. :0.0000
1st Qu.:0.000000 1st Qu.:0.00000 1st Qu.:0.00000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.000 1st Qu.:0.00000 1st Qu.:0.000 1st Qu.:0.00000 1st Qu.:0.00000 1st Qu.:0.00000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.000000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.2500 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.00000 1st Qu.:0.000000 1st Qu.:0.000 1st Qu.:0.2500 1st Qu.:0.0000 1st Qu.:0.00000 1st Qu.:0.00000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.00000 1st Qu.:1.0000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.000 1st Qu.:1.0000 1st Qu.:0.0000
Median :0.000000 Median :0.00000 Median :0.00000 Median :0.1667 Median :0.0000 Median :0.0000 Median :0.000 Median :0.00000 Median :0.000 Median :0.00000 Median :0.00000 Median :0.00000 Median :0.3000 Median :0.2500 Median :0.000000 Median :0.5000 Median :0.5000 Median :0.7500 Median :0.2500 Median :0.0000 Median :0.3333 Median :0.0000 Median :0.00000 Median :0.000000 Median :0.000 Median :0.5000 Median :0.2500 Median :0.00000 Median :0.00000 Median :0.0000 Median :1.0000 Median :0.0000 Median :0.0000 Median :0.00000 Median :1.0000 Median :0.0000 Median :0.0000 Median :0.0000 Median :0.000 Median :1.0000 Median :0.0000
Mean :0.002991 Mean :0.03174 Mean :0.07874 Mean :0.4137 Mean :0.2191 Mean :0.1678 Mean :0.086 Mean :0.17767 Mean :0.218 Mean :0.02404 Mean :0.06415 Mean :0.04167 Mean :0.2606 Mean :0.2086 Mean :0.005342 Mean :0.5741 Mean :0.4259 Mean :0.6624 Mean :0.3376 Mean :0.2485 Mean :0.4201 Mean :0.2975 Mean :0.02825 Mean :0.005698 Mean :0.251 Mean :0.4471 Mean :0.2292 Mean :0.07265 Mean :0.02137 Mean :0.2233 Mean :0.5203 Mean :0.1175 Mean :0.1175 Mean :0.03139 Mean :0.9686 Mean :0.0113 Mean :0.2096 Mean :0.4271 Mean :0.352 Mean :0.8815 Mean :0.1185
3rd Qu.:0.000000 3rd Qu.:0.00000 3rd Qu.:0.00000 3rd Qu.:1.0000 3rd Qu.:0.3333 3rd Qu.:0.0000 3rd Qu.:0.000 3rd Qu.:0.09375 3rd Qu.:0.250 3rd Qu.:0.00000 3rd Qu.:0.00000 3rd Qu.:0.00000 3rd Qu.:0.5000 3rd Qu.:0.4000 3rd Qu.:0.000000 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:0.7500 3rd Qu.:0.5000 3rd Qu.:0.7500 3rd Qu.:0.5000 3rd Qu.:0.00000 3rd Qu.:0.000000 3rd Qu.:0.575 3rd Qu.:0.6000 3rd Qu.:0.5000 3rd Qu.:0.00000 3rd Qu.:0.00000 3rd Qu.:0.0000 3rd Qu.:1.0000 3rd Qu.:0.0000 3rd Qu.:0.0000 3rd Qu.:0.00000 3rd Qu.:1.0000 3rd Qu.:0.0000 3rd Qu.:0.5000 3rd Qu.:1.0000 3rd Qu.:0.500 3rd Qu.:1.0000 3rd Qu.:0.0000
Max. :0.500000 Max. :1.00000 Max. :1.00000 Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :1.000 Max. :1.00000 Max. :1.000 Max. :1.00000 Max. :0.75000 Max. :1.00000 Max. :0.7500 Max. :0.5000 Max. :1.000000 Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :0.25000 Max. :1.000000 Max. :1.000 Max. :1.0000 Max. :0.7500 Max. :1.00000 Max. :1.00000 Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :1.00000 Max. :1.0000 Max. :0.5000 Max. :1.0000 Max. :1.0000 Max. :1.000 Max. :1.0000 Max. :1.0000
NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA’s :3 NA’s :3 NA’s :17 NA’s :17 NA’s :17 NA’s :17 NA’s :21 NA’s :21
## # A tibble: 10 x 2
##    Trait                                    w_tot
##    <fct>                                    <dbl>
##  1 Maximum size (mm)                            1
##  2 Feeding method                               1
##  3 Food size                                    1
##  4 "Preferred substrate position (adults) "     1
##  5 Living habit                                 1
##  6 Adult movement capacity (daily)              1
##  7 Bioturbation                                 1
##  8 "Sexual differenciation "                    1
##  9 Development mode                             1
## 10 Reproduction frequency                       1

iii) Missing data imputation

For the reproduction frequency, development mode, and sexual differentiation, data were missing for 9% (21 species), 7% (17), and 1% (3) of the species respectively. For these traits, we imputed missing values using nearest neighbour imputation relying on Gower dissimilarity that accommodates missing data. The imputation approach used in this article was originally developped for continuous traits (which is applicable here with our fuzzy-coded matrix) and was found to be particurlarly suitable when the trait distribution was unbalanced (Taugourdeau et al. 2014). Missing traits were imputed based on the median value of the functionally closest species for which the trait was known as well as those falling within a threshold dissimilarity of \(0.01\times d_{gower}\), with \(d_{gower}\) the dissimilarity between this closest species and the species to be inferred. This procedure gave similar results to imputation based on the 5 nearest neighbours using the kNN function. The species used to infer each missing data were then verified by experts of benthic taxonomy to ensure the ecological soundness of this imputation procedure.

Gower’s dissimilarity was computed on the relative scores (and not the raw fuzzy-coded scores) to account for the relative strengths of affinities within each traits (an affinity of 2 for a modality accompanied with an affinity of 2 for another modality of the same trait - meaning no preference between the two modalities of the trait - did not have the same meaning in our coding scheme that the same affinity of 2 accompanied with an affinity of 3 for the other modality - meaning a slight preference of the species for the second modality). Using the raw fuzzy-coded scores would not have accounted for this difference.

  • Imputation procedure
# Compute the initial functional dissimilarity among species
dist_sp <- as.matrix(gowdis(x = trait_prop, w = mod_w_tot))

# Format the fuzzy-coded matrix
trait_fuz_init <- trait_fuz %>%
  rownames_to_column("species") %>%
  gather(Modality, value, -species) %>%
  left_join(.,trait_list,by=c("Modality"="Accronyme")) %>%
  select(-Modalities,-Type)

# Initiliaze the loop with the dataframe to be filled (_estim to remember that some trait data were estimated)
trait_fuz_estim <- trait_fuz_init

# We identify the species and traits for which data are missing (note that in this code, when data are missing for one modality, they are missing for all other modalities of the trait) 
trait_fuz_na <- trait_fuz_init %>%
  filter(is.na(value)) %>%
  distinct(species,Trait,value) %>%
  arrange(species)

# Identify the target species
target <- unique(trait_fuz_na$species)

# Object to record all modification that are done
logfile <- NULL

# Missing traits were imputed based on the median value of the functionally closest species for which the trait was known as well as those falling within a threshold of 1% of the dissimilarity between this closest species and the species to be inferred.
thresh <- 0.01

for ( i in target){
  
  # 1) For each target species, we identify which traits have missing data
  #------------------------------------------------------------------------

  trait_miss <- trait_fuz_na %>% filter(species==i) %>% distinct(Trait) %>% pull()
  
  # 2) We identify the functionnaly similar species
  #------------------------------------------------

  # Sort the distances of the target species with all other species
  sim_sp <- dist_sp[i,]
  
  # Remove the distance of the target species with itself
  sim_sp <- sim_sp[names(sim_sp)!=i]
  
  # Identify the closest species
  prox_sp <- sort(sim_sp)[1]
  
  # Identify the other species whose distance with the target species is within 1% of the distance of the closet species
  dist_sp_to_use <- sim_sp[ sim_sp <= prox_sp + prox_sp * thresh ]
  
  # Retrieve the names of the species whose value will be used
  sp_to_use <- names(dist_sp_to_use)
  
  # Do we have for these species at least one for which we have data for the missing traits 
  nodata <- trait_fuz_init %>%
    filter(species %in% sp_to_use & Trait %in% trait_miss) %>%
    group_by(Trait) %>%
    # If not, we put TRUE to absence, if we have, we put FALSE.
    summarise(Absence = all(is.na(value))) %>%
    pull(Absence)
  
  # If we have don't have any data for this trait for the functionnaly closest species and those around (within the threhsold), we print the table of the missing data for these species...
  if (any(nodata)){
    print("###########################################")
    print(paste("Missing data for the closest species that should have been used to estimate the missing trait data of species", i))
    print("------------------------------------------")
    print(trait_fuz_init %>% filter(species %in% sp_to_use & Trait %in% trait_miss) %>% select(-Trait) %>% spread(Modality, value))
  }
  
  # ... and we continue the procedure by looking at the following closest species (and those within 1% of its distance with the target species), and so on until we have data...
  while( any(nodata) ) {
    
    # Species already looked at
    n <- length(sp_to_use)
    
    # Identify the following closest species
    prox_sp <- sort(sim_sp)[1+n]
    
    # Identify the species within 1% of the distance of this "closest" species with the target species
    dist_sp_to_use <- sim_sp[ sim_sp <= prox_sp + prox_sp * thresh ]
    
    sp_to_use <- names(dist_sp_to_use)
    
    # Check data availability
    nodata <- trait_fuz_init %>%
      filter(species %in% sp_to_use & Trait %in% trait_miss) %>%
      group_by(Trait) %>%
      summarise(Absence = all(is.na(value))) %>%
      pull(Absence)
  }
  
  # 3) Get the median of their values (remove NA's if they are any)
  #--------------------------------------------------------------
  
  estim <- trait_fuz_init %>%
    filter(species %in% sp_to_use & Trait %in% trait_miss) %>%
    group_by(Trait, Modality) %>%
    summarise(estimate = median(value,na.rm=T))
      
  # 4) Replace the missing values (with the modalities to ensure everything is in the right order)
  #----------------------------------------------------------------------------------------
  trait_fuz_estim <- trait_fuz_estim %>%
    mutate(value=replace(value,species==i & Trait %in% estim$Trait, estim$estimate), Modality=replace(Modality,species==i & Trait %in% estim$Trait, estim$Modality))
    
 # 5) We save a record of the species used for the imputation, their distance to the target species and their trait values 
  #----------------------------------------------------------------------------------------  
  logfile <- trait_fuz_init %>%
    filter(species %in% sp_to_use & Trait %in% trait_miss) %>% 
    select(-Trait) %>%
    spread(Modality, value) %>%
    mutate(target=i, DISTANCE = dist_sp_to_use, Missing=paste(as.list(as.character(trait_miss)), collapse=" / ")) %>%
    select(target,species,DISTANCE,everything()) %>%
    bind_rows(logfile,.)
}
## [1] "###########################################"
## [1] "Missing data for the closest species that should have been used to estimate the missing trait data of species Asclerocheilus intermedius"
## [1] "------------------------------------------"
##                 species Dev_asex Dev_direct Dev_lecitho Dev_plankto
## 1 Sclerocheilus minutus       NA         NA          NA          NA
##   Iteroparous Semelparous
## 1          NA          NA
## [1] "###########################################"
## [1] "Missing data for the closest species that should have been used to estimate the missing trait data of species Goniadella spp."
## [1] "------------------------------------------"
##                species Dev_asex Dev_direct Dev_lecitho Dev_plankto
## 1  Labioleanira yhleni       NA         NA          NA          NA
## 2 Neoleanira tetragona       NA         NA          NA          NA
## 3    Pelogenia arenosa       NA         NA          NA          NA
##   Gonochoric Hermaphrodite Iteroparous Semelparous
## 1          4             0           4           0
## 2          4             0           4           0
## 3          4             0           4           0
## [1] "###########################################"
## [1] "Missing data for the closest species that should have been used to estimate the missing trait data of species Neoleanira tetragona"
## [1] "------------------------------------------"
##             species Dev_asex Dev_direct Dev_lecitho Dev_plankto
## 1 Pelogenia arenosa       NA         NA          NA          NA
## [1] "###########################################"
## [1] "Missing data for the closest species that should have been used to estimate the missing trait data of species Nereiphylla rubiginosa"
## [1] "------------------------------------------"
##                species Iteroparous Semelparous
## 1 Notophyllum foliosum          NA          NA
## [1] "###########################################"
## [1] "Missing data for the closest species that should have been used to estimate the missing trait data of species Notophyllum foliosum"
## [1] "------------------------------------------"
##                  species Iteroparous Semelparous
## 1 Nereiphylla rubiginosa          NA          NA
## [1] "###########################################"
## [1] "Missing data for the closest species that should have been used to estimate the missing trait data of species Pelogenia arenosa"
## [1] "------------------------------------------"
##                species Dev_asex Dev_direct Dev_lecitho Dev_plankto
## 1 Neoleanira tetragona       NA         NA          NA          NA
## [1] "###########################################"
## [1] "Missing data for the closest species that should have been used to estimate the missing trait data of species Pilargis verrucosa"
## [1] "------------------------------------------"
##                species Dev_asex Dev_direct Dev_lecitho Dev_plankto
## 1 Neoleanira tetragona       NA         NA          NA          NA
## 2    Pelogenia arenosa       NA         NA          NA          NA
##   Iteroparous Semelparous
## 1           4           0
## 2           4           0
## [1] "###########################################"
## [1] "Missing data for the closest species that should have been used to estimate the missing trait data of species Polycirrus spp."
## [1] "------------------------------------------"
##                  species Iteroparous Semelparous
## 1 Amphitritides gracilis          NA          NA
## [1] "###########################################"
## [1] "Missing data for the closest species that should have been used to estimate the missing trait data of species Scalibregma celticum"
## [1] "------------------------------------------"
##                      species Dev_asex Dev_direct Dev_lecitho Dev_plankto
## 1 Asclerocheilus intermedius       NA         NA          NA          NA
## 2      Sclerocheilus minutus       NA         NA          NA          NA
## [1] "###########################################"
## [1] "Missing data for the closest species that should have been used to estimate the missing trait data of species Sclerocheilus minutus"
## [1] "------------------------------------------"
##                      species Dev_asex Dev_direct Dev_lecitho Dev_plankto
## 1 Asclerocheilus intermedius       NA         NA          NA          NA
##   Iteroparous Semelparous
## 1          NA          NA
## [1] "###########################################"
## [1] "Missing data for the closest species that should have been used to estimate the missing trait data of species Sigalion mathildae"
## [1] "------------------------------------------"
##                species Dev_asex Dev_direct Dev_lecitho Dev_plankto
## 1 Neoleanira tetragona       NA         NA          NA          NA
## 2    Pelogenia arenosa       NA         NA          NA          NA
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.02679 0.03150 0.03699 0.04677 0.13439

The mean distance of the species used to estimate missing data is 0.04. At most, the imputation was made with a species with a dissimilarity of 0.13 with the species to impute. All these functional distances are fairly small compared to those observed across all the dataset (see below). In addition, they are in line with the \(d_{gower} \leq 0.05\) criteria used by Taugourdeau et al. (2014)

  • Comparison with kNN
## [1] "Number of different imputations between our approach and kNN: 24/116"
Species Modality Article_imputation kNN_imputation
Labioleanira yhleni Dev_plankto 0 4
Streblosoma bairdi Dev_lecitho 0 4
Streblosoma bairdi Iteroparous 0 4
Notomastus latericeus Semelparous 0 4
Sphaerodoridae spp. Semelparous 0 4
Arabella iricolor Dev_plankto 0 2
Amphitritides gracilis Dev_lecitho 2 4
Claparedepelogenia inclusa Dev_lecitho 0 2
Drilonereis filum Dev_lecitho 0 2
Goniadella spp. Dev_lecitho 0 2
Lacydonia miranda Dev_lecitho 2 4
Pilargis verrucosa Dev_lecitho 0 2
Amphitritides gracilis Dev_direct 2 0
Arabella iricolor Dev_direct 2 0
Lacydonia miranda Dev_direct 2 0
Claparedepelogenia inclusa Dev_plankto 4 2
Drilonereis filum Dev_plankto 4 2
Goniadella spp. Dev_plankto 4 2
Pilargis verrucosa Dev_plankto 4 2
Streblosoma bairdi Dev_direct 4 0
Labioleanira yhleni Dev_lecitho 4 0
Notomastus latericeus Iteroparous 4 0
Sphaerodoridae spp. Iteroparous 4 0
Streblosoma bairdi Semelparous 4 0
  • Format the imputed matrix
  • Compare the trait-dissimilarity matrix computed from the initial and imputed trait matrix

iv) Assemblage‐by‐trait matrix

There are two ways to compute the assemblage-by-trait matrix containing the abundance of each modality in each site (also referred to as absolute mean value of modality, Fig.4 of Beauchard et al. 2017) :

  1. One can first compute the \(CWM\) (Community Weighted Mean) matrix through the functcomp function of the FD package for instance. The CWM matrix represents the abundance-weighted mean affinity of a site to the trait modalities, computed as \(Mean~affinity = \frac{\sum Affinity_{[i]} \times Abondance_{[i]}}{Abondance_{Tot}}\). The latter can then be multiplied by the total abundance in each site to have the abundance of the modalities.

  2. A more straightforward approach would be to do the matrix product of the site-by-species matrix with the species-by-trait matrix.

Both are be equivalent when there are no NA’s, which is the case given that we imputed all missing values. However, note that both approach would differ in the presence of NA’s. In the first approach, the CWM exclude NA’s in its computation. When multiplied back with the total abundance of each site (and if including in this total abundance the species for which trait data are missing), all traits will have the same total abundance (and missing data would then be considered as mean values). In contrast, with the second approach, the total abundance of the traits (and therefore their potential weights in some multivariate analyses) would vary according to their number of missing values (NA’s are just discarded). Note also that the functcomp function consider modalities with only 0 and 1 as binary by defaults. To treat these variables as numeric (and therefore take them into account in the CWM), we need to use the bin.num argument.

## [1] TRUE
## [1] TRUE

v) Preliminary analyses (Not included in the article)

  • Correlation between modalities

The traits and modalities chosen does not appear redundant

The most correlated modalities are :

  • Inevitable within-trait trade-offs :
    • Micro/macrophageous
    • Infauna/Epifauna
    • Hermaphrodite/Gonochoric
    • Direct/planktonic developtment
  • Inevitable across-trait associations :
    • Parasitic/attached
  • Predator & scavenger with macrophagous and epibenthic lifestyle
  • Association between modalities and taxonomy (at the Family level)

First, the number of species of each Family expressing each modality of a trait was computed to get a contingency table \(Families~\times~Modalities\)

Then, the contribution of each cell of this contingency table to the \(\chi^2\) was computed as:

\(Contribution~\chi^2=\frac{(Observed-Expected)^2}{Expected}\) This was done for each trait separately

# Compute presence/absence of family/modality association
trait_fam <- trait_mat %>% 
  select(-Kingdom,-Phylum,-Class,-Order,-Genus,-Species) %>%
  gather(trait,value,-Family,-species) %>%
  # Transform into presence/absence and replace NA's by absences
  mutate(value=if_else(value > 0, 1, 0,0)) %>%
  # Compute the number of species of each family associated to each modality
  group_by(Family,trait) %>%
  summarise(value = sum(value)) %>%
  ungroup() %>%
  # Build the FamilyxModality contingency table
  spread(trait,value) %>%
  # Reorder and remove the Bioturb_R modality
  select(Family,trait_list$Accronyme) %>%
  as.data.frame()

# Set the rownames
trait_fam_c <- trait_fam %>%
  set_rownames(.$Family) %>%
  select(-Family)

# Contributions are computed for each trait separately
khi_contrib <- NULL

for(i in 1:n_trait){
  # Select the modalities of the trait
  tmp <- trait_fam_c[,seq(from=n_mod[i]+1, to=n_mod[i+1], by=1)]
  
  # Compute expected values
  khi <- chisq.test(tmp)

  # Compute the contribution to the khi^2
  tmp <- ((khi$observed - khi$expected)^2) / khi$expected
  
  # Save
  khi_contrib <- cbind(khi_contrib,tmp)
}

khi_contrib <- khi_contrib %>%
  as.data.frame() %>%
  mutate(Family=rownames(.)) %>%
  gather(Modality,value,-Family) %>%
  left_join(.,trait_list,by=c("Modality"="Accronyme"))

# Plot
p <- ggplot(khi_contrib,aes(x=Modality,y=Family,fill=value))
p <- p + geom_tile(col="white",alpha=0.7)
p <- p + facet_grid(Family~Trait,scales="free",space="free",switch="y",labeller=label_wrap_gen(width = 10, multi_line = TRUE))
p <- p + scale_fill_gradient(low = "white", high = "black",na.value = "white",name=expression(Contribution~chi^2))
p <- p + theme_light() + xlab("")
p <- p + theme(text=element_text(size=24),axis.text.x=element_text(angle=90,hjust = 1),strip.text.x=element_text(angle=90),strip.text.y=element_text(angle=45),strip.text=element_text(size=26,face="bold"),axis.text.y=element_blank(),axis.ticks.y=element_blank())

plot(p)


Metadata (environment & coordinates)

# Environmental data
envir <- depth_fetch %>%
  select(-latitude,-longitude) %>%
  left_join(.,previmer) %>%
  left_join(.,granulo) %>%
  # Remove the site of Keraliou with incomplete data
  filter(site!="Rade de Brest - Keraliou") %>%
  # Clean site names and order the sites and habitats as in the community data
  mutate(site=recode(site,"Saint-Brieuc"="Baie de Saint-Brieuc")) %>%
  mutate(site=ifelse(site=="Baie de Saint-Brieuc" & theme=="Bancs de Maërl","Arcouest",site)) %>%
  # Translate the habitats
  mutate(theme=recode(theme, "Bancs de Maërl" = "Subtidal maerl beds", "Subtidal Meuble" = "Subtidal bare sediments", "Herbiers Intertidaux" = "Intertidal seagrass beds", "Intertidal Meuble" = "Intertidal bare sediments")) %>%
  left_join(id_pol,.) %>%
  rename(Mud = mud)

# Coordinates
coord <- coord %>%
  # Change the name of the sites as in the community data
  mutate(site = case_when(
    str_detect(point_name,"Rozegat") ~ "Rade de Brest - Rozegat",
    str_detect(point_name,"Keraliou") ~ "Rade de Brest - Keraliou",
    TRUE ~ as.character(site_name)
  )) %>%
  # Remove the site of Keraliou with incomplete data
  filter(site!="Rade de Brest - Keraliou") %>%
  # Rename the sites
  mutate(site=recode(site,"Saint-Brieuc"="Baie de Saint-Brieuc")) %>%
  mutate(site=ifelse(site=="Baie de Saint-Brieuc" & theme=="Bancs de Maërl","Arcouest",site)) %>%
  # Select only the central point of the sites as well as the only point in Pierre Noire
  filter(str_detect(point_name,"2") | str_detect(point_name,".*B$") | point_name=="PN") %>%
  select(theme,site,longitude,latitude) %>% 
  # Translate the habitats
  mutate(theme=recode(theme, "Bancs de Maërl" = "Subtidal maerl beds", "Subtidal Meuble" = "Subtidal bare sediments", "Herbiers Intertidaux" = "Intertidal seagrass beds", "Intertidal Meuble" = "Intertidal bare sediments"))

Summary of the data object

  • Metadata

    • id_tot : sample metadata for the whole community data
    • id_pol : sample metadata for the polychaetes data
  • Community matrix (site-by-species matrix)

    • com_tot : whole community matrix (Dimensions : 148, 712)
    • com_pol : community matrix with only polychaetes (Dimensions : 148, 234)
    • com_pol_inter : community matrix with only the polychaetes for the intertidal (Intertidal bare sediments and Intertidal seagrass beds;Dimensions : 79, 134)
    • com_pol_sub : community matrix with only the polychaetes for the subtidal (Subtidal bare sediments and Subtidal maerl beds;Dimensions : 69, 221)
  • Species-by-trait matrix

    • trait_list : the list of the traits and associated modalities

    • trait_fuz : trait matrix fuzzy coded (Dimensions : 234, 41)
    • trait_prop : trait matrix transformed for all modalities of each trait to sum to 1 (Dimensions : 234, 41)
    • trait_fuz_estim : trait matrix fuzzy coded with imputed missing data
    • trait_prop_estim : the imputed trait matrix transformed for all modalities of each trait to sum to 1 (Dimensions : 234, 41)
    • mod_w_tot : vector containing the weights to be given to each modality so that all traits have equal weights in the analyses at the end, irrespective of their number of associated modalities

  • Trait abundance matrix (site-by-trait matrix)

    • com_trait_estim : Site-by-trait matrix computed from the imputed species-by-trait matrix (Dimensions : 148, 41)

Figure 1

Figure 1 was made on Inkscape with the following R outputs

Brittany

# Format the data 
#----------------

map_coord <- id_pol %>%
  distinct(theme,site) %>%
  full_join(., coord) %>%
  mutate(site=factor(site,levels=c("Baie du Mont Saint-Michel","Saint-Benoit","Saint-Malo","Saint-Briac","Saint-Cast","Baie de Saint-Brieuc","Arcouest","Sept-Iles","Lannion","Saint-Efflam","Pierre Noire","Morlaix","Callot","Sainte-Marguerite","Molene","Blancs Sablons","Iroise","Camaret","Roscanvel","Rade de Brest",
"Rade de Brest - Rozegat","Rade de Brest - Larmor","Baie de Douarnenez Nord","Plage de l'Aber","Baie de Douarnenez Sud","Audierne","Glenan","Concarneau","Trevignon","Gavres","Erdeven","Lorient Etel","Plouharnel","Quiberon","Belle-ile","Meaban","Kerjouanno","Arradon","Vilaine Large Sud","Vilaine Large Nord","Vilaine Cote","Damgan"),ordered=T)) %>%
  mutate(site_number=as.numeric(site)) %>%
  as.data.frame(.)

# Map
#------

# Retrieve a map layer for France
france_data <- map_data("france")
names(france_data) <- c("X","Y","PID","POS","region","subregion")

# Zoom on Brittany
xlim <- c(-5.5,-1)
ylim <- c(46.5,49.5)
francemap <- clipPolys(france_data, xlim=xlim,ylim=ylim, keepExtra=TRUE)

p <- ggplot(data=map_coord,mapping=aes(y=latitude,x=longitude,fill=theme)) 
p <- p + coord_map(xlim=xlim,ylim=ylim)
p <- p + geom_polygon(data=france_data,aes(X,Y,group=PID),fill="gray85",inherit.aes=F,show.legend=F) 
p <- p + geom_path(data=france_data,aes(X,Y,group=PID),inherit.aes=F,show.legend=F, col="gray85") 
p <- p + ylab("")+ xlab("")
p <- p +scale_y_continuous(breaks=c(46.5,47,47.5,48,48.5,49,49.5,50), labels = c("46.5°N","47ºN","47.5°N", "48ºN","48.5°N", "49ºN","49.5°N","50°N")) 
p <- p + scale_x_continuous(breaks=c(-6,-5,-4,-3,-2,-1,0), labels = c("6°W","5°W","4°W","3°W","2°W","1°W","")) 
# Add the points
p <- p + geom_point(shape=21,size=4,alpha=0.5,col="black")
# Contour without transparency
p <- p + geom_point(aes(y=latitude,x=longitude),inherit.aes=F,shape=21,size=4,fill=alpha("white",0),col="black")
# Color for the habitats
p <- p + scale_fill_manual(values=c("#440154FF","#31688EFF","#35B779FF","#f1a340"),name="Habitat")
p <- p + theme_map() + theme(axis.ticks=element_line(colour="gray45"),axis.text=element_text(colour="gray45"),panel.background = element_rect(fill="aliceblue"), legend.position="none")
p <- p + guides(fill=guide_legend(override.aes=list(size=15,shape=22)))

# Add the scale
map <- p + scalebar(data=map_coord,dist=50,dist_unit="km",st.bottom=F,st.size=4,transform=TRUE,model="WGS84",anchor=c(x=-4,y=46.6))

plot(map)

Figure 2 & Figure S7

Included in the article

Not included in the article

Comparison with the whole community data


Figure 3

Here, the Euclidean distance between species was computed from the imputed trait matrix, but note that an alternative approach would have been to use the Gower dissimilarity (on the imputed or raw matrix as this dissimilarity accomodates missing values). Both approach yielded similar results.

Code to compute the Gower distance (square-rooted to retrieve euclidean properties) on the imputed matrix

We also remove assemblages with less than 5 species in order to be able to keep 5 PCoA axes (out of 31) for the calculation of the Fric, which leads to a quality of the reduced-space representation of \(0.655\).

## [1] TRUE
## FRic: Dimensionality reduction was required. The last 26 PCoA axes (out of 31 in total) were removed. 
## FRic: Quality of the reduced-space representation = 0.6553787 
## CWM: When 'x' is a distance matrix, CWM cannot be calculated.
# Format the output of dbFD
tmp_indic <- FD_tot_estim[sapply(FD_tot_estim,length,simplify=T)==nrow(id_tmp)] %>% # Keep only the object whose length is equal to the number of sites (in order to extract only the indices from the dbFD outputs)
  map_dfc(as.vector) %>%
  bind_cols(id_tmp,.) %>%
  left_join(id_pol,.) %>%
  gather(metric,value,-theme,-mini_theme,-site,-method_echant,-annee) %>%
  mutate_if(is.character,as.factor) %>%
  mutate(metric=factor(metric,levels=c("nbsp","sing.sp","FRic","FEve","FDiv","FDis","RaoQ"))) %>%
  mutate(metric=recode(metric, "RaoQ"="RaoQ_stand"))

# Compute redundancy
#--------------------

redundancy_tot_estim <- rao.diversity(com_pol,trait_prop_estim)

# Format the output : Remove the call form the list and transform into a dataframe
tmp_redund <- redundancy_tot_estim[-1] %>% 
  as.data.frame() %>%
  # Extract the rownames that identify the samples (ech for echantillon in french)
  rownames_to_column("ech") %>%
  # Compute the FRed 
  mutate(FRed = 1- (FunRao / Simpson), 
         sp_rich = as.numeric(specnumber(com_pol)), 
         abond_tot = as.numeric(rowSums(com_pol))) %>%
  # Gather
  gather(metric,value, -ech) %>% 
  # Recreate the id
  mutate(mini_theme=sapply(str_split(.$ech,"_"),'[[',1,simplify=T), site=sapply(str_split(.$ech,"_"),'[[',2,simplify=T), annee=sapply(str_split(.$ech,"_"),'[[',3,simplify=T)) %>%
  # Join with metadata to order the output
  left_join(id_pol,.) %>%
  select(-ech) %>%
  mutate(metric=factor(metric,levels=c("abond_tot","sp_rich","Simpson","FunRao","FRed","FunRedundancy"),ordered=T)) %>%
  # We remove the FunRedundancy compute as "Simpson - Rao" to keep only the one we want : 1 - Rao/Simpson
  filter(metric != "FunRedundancy") %>%
  droplevels() %>%
  mutate(metric=recode(metric, "FunRao"="RaoQ_syn")) %>%
  mutate_if(is.character,as.factor)  

Included in the article

Data summary
Name Piped data
Number of rows 296
Number of columns 7
_______________________
Column type frequency:
numeric 1
________________________
Group variables theme, metric

Variable type: numeric

skim_variable theme metric n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
value Subtidal maerl beds Richness 0 1.00 52.48 10.92 32 45.00 52.0 59.50 73 ▂▅▇▅▃
value Subtidal maerl beds Richness_assemblages_with_more_than_5_species 0 1.00 52.48 10.92 32 45.00 52.0 59.50 73 ▂▅▇▅▃
value Subtidal bare sediments Richness 0 1.00 29.40 14.42 6 19.25 26.5 38.00 68 ▃▇▅▂▁
value Subtidal bare sediments Richness_assemblages_with_more_than_5_species 0 1.00 29.40 14.42 6 19.25 26.5 38.00 68 ▃▇▅▂▁
value Intertidal seagrass beds Richness 0 1.00 24.52 8.96 10 19.00 22.0 29.00 50 ▅▇▅▂▁
value Intertidal seagrass beds Richness_assemblages_with_more_than_5_species 0 1.00 24.52 8.96 10 19.00 22.0 29.00 50 ▅▇▅▂▁
value Intertidal bare sediments Richness 0 1.00 11.60 7.51 1 4.75 10.0 15.25 29 ▇▇▆▂▃
value Intertidal bare sediments Richness_assemblages_with_more_than_5_species 14 0.73 14.68 6.39 6 10.00 14.0 18.00 29 ▇▇▂▃▂

Note that for total abundance and taxonomic richness, assemblages with <= 5 species were included, while they were removed for functional indices computation. The removal of these assemblages with <= 5 species only concerns/affects 14 assemblages of Intertidal bare sediment

# Taxnomic indices
#-----------------

taxo_div_indic <- func_indic %>%
  filter(metric %in% c("Abundance","Richness","Simpson"))

# Compute the mean value
taxo_div_indic_mean <- taxo_div_indic %>%
  group_by(theme,metric) %>%
  summarise(mu=mean(value,na.rm=T)) %>%
  # Define the limits of the segment we want to draw
  mutate(ymu = case_when(
   theme=="Subtidal maerl beds" ~ 1.4,
   theme=="Subtidal bare sediments" ~ 2.4,
   theme=="Intertidal seagrass beds" ~ 3.4,
   theme=="Intertidal bare sediments" ~ 4.4
  ))

# Plot
p_taxo <- ggplot(taxo_div_indic,aes(x=value, y = theme, fill=theme, height = ..density..)) 
p_taxo <- p_taxo +  geom_density_ridges(aes(point_color=theme),alpha = .5, scale=0.7, point_alpha = .5, jittered_points = TRUE, position = position_raincloud(adjust_vlines = TRUE, height=0.2))
p_taxo <- p_taxo + facet_wrap(~metric,scales="free_x", ncol=4)
p_taxo <- p_taxo + geom_segment(data=taxo_div_indic_mean,aes(x=mu,xend=mu,y=theme,yend=ymu),inherit.aes=F)
p_taxo <- p_taxo + geom_point(data=taxo_div_indic_mean,aes(x=mu,y=ymu, fill=theme), shape= 21, size=2.5 ,inherit.aes=F)
p_taxo <- p_taxo + scale_fill_manual(values=c("#440154FF","#31688EFF","#35B779FF","#f1a340"), guide = FALSE)
p_taxo <- p_taxo + scale_discrete_manual("point_color", values=c("#440154FF","#31688EFF","#35B779FF","#f1a340"), guide = FALSE)
p_taxo <- p_taxo + ylab("") + xlab("")  + theme_light() + theme(panel.grid.minor = element_blank(), strip.text=element_text(face="bold",colour="black"), text=element_text(size=16))


# Functional indices
#-------------------

func_div_indic <- func_indic %>%
  filter(metric %in% c("FRic","FEve","FDiv","FDis"))

# Compute the mea,
func_div_indic_mean <- func_div_indic %>%
  group_by(theme,metric) %>%
  summarise(mu=mean(value,na.rm=T)) %>%
  # Define the limits of the segment we want to draw
  mutate(ymu = case_when(
   theme=="Subtidal maerl beds" ~ 1.4,
   theme=="Subtidal bare sediments" ~ 2.4,
   theme=="Intertidal seagrass beds" ~ 3.4,
   theme=="Intertidal bare sediments" ~ 4.4
  ))

# Représentation
p_func <- ggplot(func_div_indic,aes(x=value, y = theme, fill=theme, height = ..density..)) 
p_func <- p_func + geom_density_ridges(aes(point_color=theme),alpha = .5, scale=0.7, point_alpha = .5, jittered_points = TRUE, position = position_raincloud(adjust_vlines = TRUE, height=0.2))
p_func <- p_func + facet_wrap(~metric,scales="free_x", ncol=4)
p_func <- p_func + geom_segment(data=func_div_indic_mean,aes(x=mu,xend=mu,y=theme,yend=ymu),inherit.aes=F)
p_func <- p_func + geom_point(data=func_div_indic_mean,aes(x=mu,y=ymu, fill=theme), shape= 21, size=2.5 ,inherit.aes=F)
p_func <- p_func + scale_fill_manual(values=c("#440154FF","#31688EFF","#35B779FF","#f1a340"), guide = FALSE)
p_func <- p_func + scale_discrete_manual("point_color", values=c("#440154FF","#31688EFF","#35B779FF","#f1a340"), guide = FALSE)
p_func <- p_func + ylab("") + xlab("")  + theme_light() + theme(panel.grid.minor = element_blank(), strip.text=element_text(face="bold",colour="black"), text=element_text(size=16))

grid.arrange(p_taxo,p_func, ncol=3, widths = c(1, 10, 1), heights=c(1,1), layout_matrix = rbind(c(NA, 1, NA), c(2, 2, 2))) # export size : 1000*1200

Not included in the article

  • Distribution of all indices explored during this study


Figure 4

Computation

Null model

# Initialisation
#---------------

set.seed(1)
null_model<- "trialswap"
nperm <- 1000

mat_trait <- trait_prop_estim

res <- NULL

# Loop for each permutation
#--------------------------

for(i in seq(nperm)){
  
    # Permute species separately within the intertidal and within the subtidal
    com_inter_i <- randomizeMatrix(com_pol_inter, null.model = null_model, iterations = 100000)
    com_sub_i <- randomizeMatrix(com_pol_sub, null.model = null_model, iterations = 100000)
    
    # Compute Rao's Q
    fd_inter_i <- rao.diversity(com_inter_i,mat_trait)
    fd_sub_i <- rao.diversity(com_sub_i,mat_trait)
    
    # Save the results
    res_inter_i <- data.frame(Permutation = i, Samples = rownames(com_inter_i), RaoQ = fd_inter_i$FunRao)
    res_sub_i <- data.frame(Permutation = i, Samples = rownames(com_sub_i), RaoQ = fd_sub_i$FunRao)
        
    res <- bind_rows(res,res_inter_i,res_sub_i)
}


# Format the output
#------------------

null_div_swap <- res %>% 
  # On recréer les colonnes d'identification
  mutate(mini_theme=sapply(str_split(.$Samples,"_"),'[[',1,simplify=T), site=sapply(str_split(.$Samples,"_"),'[[',2,simplify=T), annee=sapply(str_split(.$Samples,"_"),'[[',3,simplify=T)) %>%
  mutate_if(is.character,as.factor)

# Compute the mean and sd across the randomized communities for each site and year (annee in French), and habitat (mini_theme here) 
#-----------------------------------------------------------------------------------

null_div_ses_swap <- null_div_swap %>%
  group_by(mini_theme,site,annee) %>%
  summarise(n = n(), RaoQ_mu = mean(RaoQ), RaoQ_sd = sd(RaoQ))

Figure of the article

Figure 4b

# Merge SES with the site coordinates
map_ses <- coord %>%
  right_join(., ses) %>%
  mutate_if(is.character,as.factor) %>%
  mutate(mini_theme = factor(mini_theme,levels=c("IM","HI","SM","MA"), ordered = T))

# Map
#----


# Retrieve a map layer for France
france_data <- map_data("france")
names(france_data) <- c("X","Y","PID","POS","region","subregion")

# Zoom on Brittany
xlim <- c(-5.5,-1)
ylim <- c(46.5,49.5)
francemap <- clipPolys(france_data, xlim=xlim,ylim=ylim, keepExtra=TRUE)

# Plot the map
p <- ggplot(data=map_ses,mapping=aes(y=latitude,x=longitude,fill=value)) 
p <- p + facet_wrap(mini_theme~annee, ncol=6, labeller = labeller(mini_theme = as_labeller(c( "MA" = "Subtidal maerl beds", "SM" = "Subtidal bare sediments", "HI" = "Intertidal seagrass beds", "IM" = "Intertidal bare sediments"), label_wrap_gen(20,multi_line=T))))
p <- p + coord_map(xlim=xlim,ylim=ylim)
p <- p + geom_polygon(data=france_data,aes(X,Y,group=PID),fill="gray90",col=NA,inherit.aes=F,show.legend=F) 
p <- p + ylab("")+ xlab("")
p <- p +scale_y_continuous(breaks=c(46.5,47,47.5,48,48.5,49,49.5,50), labels = c("46.5°N","47ºN","47.5°N", "48ºN","48.5°N", "49ºN","49.5°N","50°N")) 
p <- p + scale_x_continuous(breaks=c(-6,-5,-4,-3,-2,-1,0), labels = c("6°W","5°W","4°W","3°W","2°W","1°W","")) 
# Add points
p <- p + geom_point(shape=21,size=4,alpha=1,stroke = 1,col="black")
# Define the colour scale
p <- p + scale_fill_gradient2(limits=c(-2,2), low = "#ca0020", mid = "#f7f7f7", high = "#0571b0", name=expression('SES'['RaoQ']))
p <- p + theme_light() + theme(text = element_text(size=16),strip.text=element_text(face="bold",colour="black"),strip.text.y=element_text(angle=0))

plot(p)


Figure 5

Included in the article

Not included in the article

Figure 6 and associated text

Coinertia of Figure 6

## [1] "The RV coefficient between the two ordinations is 0.62"
# Plot of Figure 6
#-----------------

# Extract the coordinates :

# (1) of the sites in the trait space
pos_site_trait <- coin_glob$mY %>%
rownames_to_column("location") %>%
separate(location, c("mini_theme", "site", "annee"), sep = "_") %>%
select(mini_theme, site, annee, contains("S1"), contains("S2")) %>%
rename("x" = "NorS1", "y" = "NorS2") %>%
mutate(component = "Trait")

# (2) of the sites in the taxonomic space
pos_site_taxo <- coin_glob$mX %>%
rownames_to_column("location") %>%
separate(location, c("mini_theme", "site", "annee"), sep = "_")   %>%
select(mini_theme, site, annee, contains("S1"), contains("S2")) %>%
rename("x" = "NorS1", "y" = "NorS2") %>%
mutate(component = "Taxonomy")

# (3) of the trait modalities
pos_mod_trait <- coin_glob$l1

# Format the site coordinates
coin_plot_data <- bind_rows(pos_site_trait, pos_site_taxo) %>%
  full_join(id_pol, .) %>%
  rename(habitat = theme) %>%
  unite(sample, mini_theme, site, annee, remove = F) %>%
  mutate(habitat = fct_rev(habitat))

# COmpute the group centroids
cent_group <- coin_plot_data %>%
  group_by(habitat,component) %>%
  summarise(x=mean(x),y=mean(y)) %>%
  unite(label, habitat, component, sep=" | ",remove=FALSE)


coin_plot_data_background <- coin_plot_data %>%
  rename(habitat_2 = habitat)


p <- ggplot(coin_plot_data, aes(x = x, y = y, col = habitat, shape = component))
p <- p + geom_point(data=coin_plot_data_background, aes(x=x,y=y, col=habitat_2), shape=16,inherit.aes = FALSE, alpha = 0.4, size=0.8)
p <- p + stat_chull(aes(fill=habitat,linetype=component), alpha=0.15, geom = "polygon")
p <- p + geom_point(data=cent_group, aes(x = x, y = y,fill=habitat),col="black", size=3, alpha= 0.7)
p <- p + geom_path(aes(group = habitat),alpha=0.3)
p <- p + facet_wrap(~habitat,ncol=2)
p <- p + scale_shape_manual(values = c(21, 24), name = "Centroid")
p <- p + scale_colour_manual( values = rev(c("#440154FF", "#31688EFF", "#35B779FF", "#f1a340")), guide = FALSE)
p <- p + scale_fill_manual( values = rev(c("#440154FF", "#31688EFF", "#35B779FF", "#f1a340")), guide = TRUE, name="Habitat")
p <- p + scale_linetype_manual(values=c(1,2), name="Convex hull")
p <- p + guides(fill = guide_legend(override.aes = list(shape=NA, alpha=0.7)),linetype = guide_legend(override.aes = list(colour="black"), nrow=1), shape = guide_legend(override.aes = list(fill="grey")))
p_coin <- p + theme_bw() + theme(strip.text.x = element_text(face="bold"),text = element_text(size = 16)) + xlab("Coinertia 1") + ylab("Coinertia 2")

plot(p_coin)

RV coefficients per habitat (in the text of the article)

## [1] "The RV coefficient within intertidal bare sediment is of 0.56"
## [1] "The RV coefficient within intertidal seagrass beds is of 0.85"
## [1] "The RV coefficient within subtidal bare sediment is of 0.71"
## [1] "The RV coefficient within subtidal maerl beds is of 0.54"

Table 1

Trait Modalities Abbreviations
Maximum size (mm) <2 Size_inf2
2 to 5 Size_2.5
5 to 10 Size_5.10
10 to 50 Size_10.50
50 to 100 Size_50.100
100 to 200 Size_100.200
>200 Size_sup200
Feeding method Subsurface deposit feeder SSDF
Surface deposit feeder SDF
Active suspension feeder ASF
Passive suspension feeder PSF
Grazer Grazer
Predator Pred
Scavenger Scav
Parasitic Parasitic
Food size Microphageous Microphagous
Macrophagous Macrophagous
Preferred substrate position (adults) Infaunal Infaunal
Epibenthic Epibenthic
Living habit Tube dweller Tube_dweller
Burrower Burrower
Crawler Crawler
Swimmer Swimmer
Attached Attached
Adult movement capacity (daily) None (0m) Mob_0
<10m Mob_inf10
10-100m Mob_10.100
100 - 1000m Mob_100.1000
Bioturbation None Bioturb_N
S Bioturb_S
B Bioturb_B
UC Bioturb_UC
DC Bioturb_DC
Sexual differenciation Hermaphrodite Hermaphrodite
Gonochoric Gonochoric
Development mode Asexual Dev_asex
Direct Dev_direct
Indirect - planktotrophic Dev_plankto
Indirect - lecithotrophic Dev_lecitho
Reproduction frequency Iteroparous Iteroparous
Semelparous Semelparous

Table 2

Computation

Trait space occupancy

The Euclidean distance was used to compute the trait space to be coherent with the method used to compute the functional indices for Figure 3

Here,we computed the trait space using a PCoA. This choice was originally done to compare the results presented in the article (obtained from the Euclidean distance on the imputed trait matrix) with those obtained from the square-root of Gower distance on the raw trait matrix (Gower distance is classically used to handle NAs, and was square-rooted to make it fully embeddable into Euclidean space and avoid negative eigenvalues in the PCoA;see p.297 of Legendre and Legendre 2012). However, given that we used the Euclidean distance, note that our approach is fairly equivalent to a PCA computed on the imputed trait matrix. The alternative approach with the Gower distance yielded similar results

Note that there is a typo in the legend of Table 2, as it says “trait space was measured based on the first 6 axes of the PCA of the species‐by‐trait matrix” instead of “trait space was measured based on the first 6 axes of the PCoA of the species‐by‐trait matrix”

  • How many axes should we keep ? (For computation efficiency, we cannot keep all axes otherwise the computation of the volume is fairly demanding)

We kept the 6 axes that had eigenvalues superior (or close for the 6th one) to those expected from brocken sticks. These 6 axes represented 70.36 % of the intial variance of the species-by-trait matrix.

  • Computation of the volume of the regional trait space (formed by all the species found in the survey)
  • Computation of the volume of the trait space of each habitat (formed by all the species found in these habitat across the survey)
  • Computation of the mean volume of the trait spaces of the assemblages of each habitat (space formed by the species found in each assemblage, and then average per habitat)

Given that we chose to keep 6 axes to accurately depict the functional trait space, we cannot compute the convex hulls for the assemblages with less than 7 species. Therefore, we removed the assemblages with less than 7 species from this analysis. See below for a sensitivity analysis of the results of Table 2 to this choice as well as for the identity of the samples removed (most of them belonging to intertidal bare sediments).

The code used to compute the results found in Table 2 (based on 6 PCoA axes for all assemblages with more than 7 species) is :

Table of the article

Habitat Taxonomic BDtot Functional BDtot Total occupancy of regional trait space (%) Average occupancy of regional trait space (%) ± SD Total contribution to regional taxonomic richness (%) Average contribution to regional taxonomic richness (%) ± SD
Intertidal bare sediments 0.75 0.13 61.77 2.76 ± 4.02 40.17 4.96 ± 3.21
Intertidal seagrass beds 0.52 0.06 64.14 9.33 ± 6.52 47.01 10.48 ± 3.83
Subtidal bare sediments 0.60 0.06 82.27 15.67 ± 12.88 60.26 12.57 ± 6.16
Subtidal maerl beds 0.47 0.03 86.10 28.24 ± 7.93 77.78 22.43 ± 4.67

Note that there is an error in the printed version of the article: the 6th column of Table 2 should be named “Total contribution to regional taxonomic richness (%)” instead of “Total contribution to regional taxonomic richness (%) ± SD”. It is indeed the total contribution of the habitat, which is a unique value. Standard deviation was only used when measuring the individual contribution of the assemblages of each habitat (values that are found in the 7th and last column of Table 2).

Not included in the article

Sensitivity analysis of trait space occupancy

Given that we chose to keep 6 axes to accurately depict the functional trait space, we removed the assemblages with less than 7 species for which we could not compute the convex hulls.

Most of the assemblages removed belonged to intertidal bare sediments :

Assemblages removed Number of species
IM_Blancs Sablons_2007 1
IM_Audierne_2007 1
IM_Audierne_2010 2
IM_Saint-Briac_2007 3
IM_Saint-Briac_2010 3
IM_Blancs Sablons_2013 3
IM_Audierne_2013 3
IM_Baie du Mont Saint-Michel_2010 4
IM_Baie du Mont Saint-Michel_2013 4
IM_Saint-Benoit_2013 4
IM_Saint-Briac_2013 4
IM_Saint-Efflam_2010 4
IM_Blancs Sablons_2010 4
IM_Plage de l’Aber_2013 5
SM_Audierne_2010 6
IM_Erdeven_2010 6

However, to see if we missed something by removing these assemblages, we computed the volume occupied by these assemblages based on a reduced set of coordinates (using as many axes as possible, e.g. based on 4 PCoA axes for an assemblage with 5 species). Note that these estimated volumes are less precise that those reported for the assemblages with at least 7 species, as the reduced set of PCoA axes does not fully represent the functional space anymore (remember that the 6 PCoA axes already represent less than 80% of the original variance). Furthermore, this only possible for assemblages with at least 3 species, otherwise we cannot compute a volume.

The results (below) suggest the results reported in Table 2 (for intertidal bare sediment especially) are robust to this choice of removing the assemblages with less than 7 species.

## [1] "There are 148 assemblages. We were able to compute the convex hulls for 144 of them. 4 returned an error and could not be computed"
Assemblages that did not worked at all Number of species
IM_Blancs Sablons_2007 1
IM_Audierne_2007 1
IM_Audierne_2010 2
IM_Saint-Briac_2013 4
Assemblages with less than 7 species Percentage occupancy of regional trait space Number of axes kept for computation
IM_Saint-Briac_2007 1.66 2
IM_Saint-Briac_2010 4.03 2
IM_Blancs Sablons_2013 7.91 2
IM_Audierne_2013 2.15 2
IM_Baie du Mont Saint-Michel_2010 10.04 3
IM_Baie du Mont Saint-Michel_2013 1.16 3
IM_Saint-Benoit_2013 3.04 3
IM_Saint-Efflam_2010 1.88 3
IM_Blancs Sablons_2010 1.15 3
IM_Plage de l’Aber_2013 0.78 4
SM_Audierne_2010 0.02 5
IM_Erdeven_2010 0.00 5
Habitat SD vol.rel.sd Mean number of axes used
Intertidal bare sediments 2.830963 3.810377 5.291667
Intertidal seagrass beds 9.328768 6.523469 6.000000
Subtidal bare sediments 15.293168 12.947109 5.976191
Subtidal maerl beds 28.238898 7.926320 6.000000

Sensitivity of BDtot to alternative transformations

  • Hellinger

  • Chord

  • Log.chord


Figure S1

Figure S1 was combined on Inkscape with the following R outputs

Map of the sites

# Map
#------

# Retrieve a map layer for France
france_data <- map_data("france")
names(france_data) <- c("X","Y","PID","POS","region","subregion")

# Zoom on Brittany
xlim <- c(-5.5,-1)
ylim <- c(46.5,49.5)
francemap <- clipPolys(france_data, xlim=xlim,ylim=ylim, keepExtra=TRUE)

p <- ggplot(data=map_coord,mapping=aes(y=latitude,x=longitude,fill=theme)) 
p <- p + coord_map(xlim=xlim,ylim=ylim)
p <- p + geom_polygon(data=france_data,aes(X,Y,group=PID),fill="gray85",inherit.aes=F,show.legend=F) 
p <- p + geom_path(data=france_data,aes(X,Y,group=PID),inherit.aes=F,show.legend=F, col="gray85") 
p <- p + ylab("")+ xlab("")
p <- p +scale_y_continuous(breaks=c(46.5,47,47.5,48,48.5,49,49.5,50), labels = c("46.5°N","47ºN","47.5°N", "48ºN","48.5°N", "49ºN","49.5°N","50°N")) 
p <- p + scale_x_continuous(breaks=c(-6,-5,-4,-3,-2,-1,0), labels = c("6°W","5°W","4°W","3°W","2°W","1°W","")) 
# Add the points
p <- p + geom_point(shape=21,size=4,alpha=0.5,col="black")
# Contour without transparency
p <- p + geom_point(aes(y=latitude,x=longitude),inherit.aes=F,shape=21,size=4,fill=alpha("white",0),col="black")
# Color for the habitats
p <- p + scale_fill_manual(values=c("#440154FF","#31688EFF","#35B779FF","#f1a340"),name="Habitat")
p <- p + theme_map() + theme(axis.ticks=element_line(colour="gray45"),axis.text=element_text(colour="gray45"),panel.background = element_rect(fill="aliceblue"), legend.position="none")
p <- p + guides(fill=guide_legend(override.aes=list(size=15,shape=22)))

p <- p +  geom_text_repel(aes(y=latitude,x=longitude,label = site_number),inherit.aes=F, size= 3, point.padding=unit(0.6, "lines"),segment.size = 0.5,show.legend=F)

# Add the scale
map <- p + scalebar(data=map_coord,dist=50,dist_unit="km",st.bottom=F,st.size=4,transform=TRUE,model="WGS84",anchor=c(x=-4,y=46.6))

plot(map)

Number of samples available per site/year

# Create a df with the missing samples
missing <- data.frame(theme="Intertidal bare sediments",site=c("Rade de Brest","Plage de l'Aber"),annee="2010",method_echant="Carotte",n_rep=0,Available=FALSE)

# Retrieve the number of samples available for the remaining sampling occasions
tmp  <- faune_sel %>%
  distinct(theme, site,point_name, annee, method_echant,replicat) %>%
  # Remove the maerl of Keraliou
  filter(site!="Rade de Brest - Keraliou") %>%
  # Rename the sites
  mutate(site=recode(site,"Saint-Brieuc"="Baie de Saint-Brieuc")) %>%
  mutate(site=ifelse(site=="Baie de Saint-Brieuc" & theme=="Bancs de Maërl","Arcouest",site)) %>%
  # Summarise the number of replicates per sampling occasion
  group_by(theme, site, annee, method_echant) %>%
  dplyr::summarise(n_rep=n()) %>%
  ungroup() %>%
  # Define the availability of data
  mutate(Available=TRUE) %>%
  # Add the missing samples
  bind_rows(.,missing) %>%
  # Rename and reorder the factors
  mutate(theme=recode(theme, "Bancs de Maërl" = "Subtidal maerl beds", "Subtidal Meuble" = "Subtidal bare sediments", "Herbiers Intertidaux" = "Intertidal seagrass beds", "Intertidal Meuble" = "Intertidal bare sediments")) %>%
  mutate(theme = factor(theme,levels=c("Subtidal maerl beds","Subtidal bare sediments","Intertidal seagrass beds","Intertidal bare sediments"), ordered=TRUE)) %>%
  mutate(site=factor(site,levels=c("Baie du Mont Saint-Michel","Saint-Benoit","Saint-Malo","Saint-Briac","Saint-Cast","Baie de Saint-Brieuc","Arcouest","Sept-Iles","Lannion","Saint-Efflam","Pierre Noire","Morlaix","Callot","Sainte-Marguerite","Molene","Blancs Sablons","Iroise","Camaret","Roscanvel","Rade de Brest",
"Rade de Brest - Rozegat","Rade de Brest - Larmor","Baie de Douarnenez Nord","Plage de l'Aber","Baie de Douarnenez Sud","Audierne","Glenan","Concarneau","Trevignon","Gavres","Erdeven","Lorient Etel","Plouharnel","Quiberon","Belle-ile","Meaban","Kerjouanno","Arradon","Vilaine Large Sud","Vilaine Large Nord","Vilaine Cote","Damgan"),ordered=T)) %>%
  mutate(site_number=as.numeric(site)) %>%
  arrange(site_number) %>%
  mutate(num_site=factor(paste(site_number,site,sep=" - "),levels=rev(unique(paste(site_number,site,sep=" - "))),ordered=TRUE))

# Plot
p <- ggplot(aes(x=annee,y=num_site,fill=theme,shape=Available),data=tmp)
p <- p + geom_point(size=12,alpha=0.5)
# Add the number of samples
p <- p + geom_text(aes(x=annee,y=num_site,label=n_rep),size = 7.5,fontface="bold",show.legend=FALSE)
p <- p + facet_wrap(~theme,ncol=4,labeller=label_wrap_gen(width = 10, multi_line = TRUE))
p <- p + theme_light()
p <- p + labs(y='')
p <- p + scale_shape_manual(values=c(0,22),labels=c("None","# available"), name="Availability")
p <- p + scale_x_discrete(name="Years",drop=FALSE)
p <- p + scale_fill_manual(values=c("#440154FF","#31688EFF","#35B779FF","#f1a340"),name="Habitat")
p <- p + theme(strip.text.x = element_text(face="bold", colour="black"),axis.text.x=element_text(angle=90), text = element_text(size = 32))
dispo <- p + guides(fill = guide_legend(override.aes = list(label=FALSE,shape=22),ncol = 2,order = 1), shape=guide_legend(override.aes=list(shape=c(22,22),fill=c("white","#f1a340")),order = 2))

plot(dispo)


Figure S3

# Transformation hellinger
com_pol_hell <- decostand(com_pol,"hellinger") 

# Calcul de l'ACP
com_pol_pca <- rda(com_pol_hell)

# On récupère les coordonnées des sites
site_scores <- scores(com_pol_pca,display="sites", scaling = 1) %>%
  as.data.frame() %>%
  rownames_to_column("location") %>%
  separate(location,c("mini_theme","site","annee"),sep="_") %>%
  right_join(id_pol,.) %>% 
  mutate(annee = as.numeric(annee)) %>%
  mutate(tidal = recode(method_echant, "Benne" = "Subtidal", "Carotte" = "Intertidal")) %>%
  mutate(tidal = factor(tidal, levels= c("Intertidal", "Subtidal"), ordered = T)) %>%
  # On rassemble les sites de la Rade et de Morlaix, ainsi qu'on modifie ceux de Saint-brieuc
  mutate(site=recode(site,"Saint-Brieuc"="Baie de Saint-Brieuc")) %>%
  mutate(site=ifelse(site=="Baie de Saint-Brieuc" & theme=="Bancs de Maërl","Arcouest",site)) %>%
  # On réordonne le facteur theme et metric pour la représentation
  mutate(theme = factor(theme,levels=c("Subtidal maerl beds","Subtidal bare sediments","Intertidal seagrass beds","Intertidal bare sediments"), ordered=TRUE)) %>%
  mutate(site=factor(site,levels=c("Baie du Mont Saint-Michel","Saint-Benoit","Saint-Malo","Saint-Briac","Saint-Cast","Baie de Saint-Brieuc","Arcouest","Sept-Iles","Lannion","Saint-Efflam","Pierre Noire","Morlaix","Callot","Sainte-Marguerite","Molene","Blancs Sablons","Iroise","Camaret","Roscanvel","Rade de Brest",
"Rade de Brest - Rozegat","Rade de Brest - Larmor","Baie de Douarnenez Nord","Plage de l'Aber","Baie de Douarnenez Sud","Audierne","Glenan","Concarneau","Trevignon","Gavres","Erdeven","Lorient Etel","Plouharnel","Quiberon","Belle-ile","Meaban","Kerjouanno","Arradon","Vilaine Large Sud","Vilaine Large Nord","Vilaine Cote","Damgan"),ordered=T)) %>%
  mutate(number=as.numeric(site))
  
# On récupère les centroide de chaque site
site_scores_cent <- site_scores %>%
  group_by(theme,site) %>%
  summarise(cent1 = mean(PC1), cent2 = mean(PC2)) %>%
  ungroup %>%
  left_join(site_scores,.)

# On récupère les centroides pour les labels
label <- site_scores_cent %>%
  group_by(theme,site,number) %>%
  summarise(cent1 = mean(PC1), cent2 = mean(PC2)) %>%
  ungroup

# Calculate the variance represented by each axis
var_axes <- round(com_pol_pca$CA$eig/sum(com_pol_pca$CA$eig)*100,2)

# Représentation
#---------------

# Sites
p <- ggplot(site_scores, aes(x=PC1, y=PC2, col=theme, fill = theme, group=site))
p <- p + geom_point(size=2,shape=21, alpha=0.7) + geom_segment(data = site_scores_cent, aes(x=cent1, y=cent2, xend=PC1, yend=PC2, color=theme))
p <- p + geom_label_repel(data= label, aes(x=cent1, y=cent2, label = number, fill = theme), fontface="bold",size=3, col="white", segment.colour = "black", box.padding = 0, inherit.aes=FALSE, show.legend = FALSE, alpha=0.8)
p <- p + coord_equal()
p <- p + scale_colour_manual(values=c("#440154FF", "#31688EFF", "#35B779FF", "#f1a340"),name="")
p <- p + scale_fill_manual(values=c("#440154FF", "#31688EFF", "#35B779FF", "#f1a340"),name="")
p <- p + xlab(paste("PCA1:",var_axes[1],"%")) + ylab(paste("PCA2:",var_axes[2],"%"))
p <- p + theme_bw() + theme(text=element_text(size=16), strip.text.x = element_text(face="bold"), legend.position = "bottom")
p <- p + guides(colour = guide_legend(override.aes = list(shape = 15, size=5, alpha = 0.7)))

p


Figure S4

Computation

Null model per traits

# Initialisation
#---------------

set.seed(1) 

null_model<- "trialswap"
nperm <- 1000

mat_trait <- trait_prop_estim

res <- NULL

# Loop on the traits
#--------------------

for(t in unique(trait_list$Trait)){
  
  # Select all the modalities of the trait
  mod <- trait_list$Accronyme[trait_list$Trait == t]
  
  tmp_trait <- mat_trait[,mod]
  
  # Loop for each permutation
  for(i in seq(nperm)){
  
    # Permute species separately within the intertidal and within the subtidal
    com_inter_i <- randomizeMatrix(com_pol_inter, null.model = null_model, iterations = 100000)
    com_sub_i <- randomizeMatrix(com_pol_sub, null.model = null_model, iterations = 100000)
    
    # Compute Rao's Q
    fd_inter_i <- rao.diversity(com_inter_i,tmp_trait)
    fd_sub_i <- rao.diversity(com_sub_i,tmp_trait)
    
    # Save the resuts
    res_inter_i <- data.frame(Permutation = i, Trait = t, Samples = rownames(com_inter_i), RaoQ = fd_inter_i$FunRao)
    res_sub_i <- data.frame(Permutation = i, Trait = t, Samples = rownames(com_sub_i), RaoQ = fd_sub_i$FunRao)
        
    res <- bind_rows(res,res_inter_i,res_sub_i)
  }
}

# Format the output
#--------------

null_div_swap <- res %>%
  # Recreate the samples id
  mutate(mini_theme=sapply(str_split(.$Samples,"_"),'[[',1,simplify=T), site=sapply(str_split(.$Samples,"_"),'[[',2,simplify=T), annee=sapply(str_split(.$Samples,"_"),'[[',3,simplify=T)) %>%
  mutate_if(is.character,as.factor)

# Compute the mean and sd across the randomized communities for each site and year (annee in French), and habitat (mini_theme here)  
#---------------------------------------------------------------------------------

null_div_ses_swap <- null_div_swap %>%
  group_by(Trait,mini_theme,site,annee) %>%
  summarise(n = n(), RaoQ_mu = mean(RaoQ), RaoQ_sd = sd(RaoQ))

Figure of the article


Session info

## R version 3.5.0 (2018-04-23)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS  10.14.6
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] fr_FR.UTF-8/fr_FR.UTF-8/fr_FR.UTF-8/C/fr_FR.UTF-8/fr_FR.UTF-8
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] ggtern_3.3.0         gridExtra_2.3        cowplot_0.9.2       
##  [4] kableExtra_1.1.0     knitr_1.24           ggsn_0.5.0          
##  [7] PBSmapping_2.72.1    ggpubr_0.1.7         ggcorrplot_0.1.3    
## [10] ggthemes_4.2.0       viridis_0.5.1        viridisLite_0.3.0   
## [13] ggalt_0.4.0          ggjoy_0.4.1          ggridges_0.5.1      
## [16] ggpointdensity_0.1.0 ggrepel_0.8.2        picante_1.8         
## [19] nlme_3.1-137         adespatial_0.3-7     SYNCSA_1.3.4        
## [22] FD_1.0-12            geometry_0.4.5       ape_5.3             
## [25] ade4_1.7-13          vegan_2.5-6          lattice_0.20-35     
## [28] permute_0.9-5        skimr_2.0.2          VIM_4.8.0           
## [31] data.table_1.12.2    colorspace_1.4-1     naniar_0.4.2        
## [34] magrittr_1.5         forcats_0.3.0        stringr_1.4.0       
## [37] dplyr_0.8.5          purrr_0.3.2          readr_1.1.1         
## [40] tidyr_1.0.0          tibble_2.1.3         ggplot2_3.3.0       
## [43] tidyverse_1.2.1     
## 
## loaded via a namespace (and not attached):
##   [1] questionr_0.7.0           proto_1.0.0              
##   [3] tidyselect_0.2.5          htmlwidgets_1.3          
##   [5] ranger_0.11.2             maptools_0.9-2           
##   [7] RNeXML_2.3.0              munsell_0.5.0            
##   [9] codetools_0.2-15          units_0.6-0              
##  [11] miniUI_0.1.1.1            withr_2.1.2              
##  [13] highr_0.7                 uuid_0.1-2               
##  [15] rstudioapi_0.9.0          bayesm_3.1-4             
##  [17] robustbase_0.93-1         vcd_1.4-4                
##  [19] Rttf2pt1_1.3.7            labeling_0.3             
##  [21] RgoogleMaps_1.4.3         repr_1.0.2               
##  [23] coda_0.19-3               LearnBayes_2.15.1        
##  [25] vctrs_0.2.0               generics_0.0.2           
##  [27] xfun_0.9                  adegenet_2.1.1           
##  [29] adephylo_1.1-11           R6_2.4.1                 
##  [31] rmdformats_0.3.5          RcppArmadillo_0.9.700.2.0
##  [33] assertthat_0.2.1          promises_1.0.1           
##  [35] scales_1.0.0              nnet_7.3-12              
##  [37] gtable_0.3.0              phylobase_0.8.6          
##  [39] ash_1.0-15                rlang_0.4.0              
##  [41] zeallot_0.1.0             splines_3.5.0            
##  [43] extrafontdb_1.0           lazyeval_0.2.2           
##  [45] acepack_1.4.1             checkmate_1.9.1          
##  [47] broom_0.5.2               yaml_2.2.0               
##  [49] reshape2_1.4.3            abind_1.4-5              
##  [51] modelr_0.1.5              backports_1.1.4          
##  [53] httpuv_1.5.2              Hmisc_4.2-0              
##  [55] extrafont_0.17            tensorA_0.36.1           
##  [57] tools_3.5.0               bookdown_0.13            
##  [59] spData_0.2.9.0            ellipsis_0.3.0           
##  [61] RColorBrewer_1.1-2        latex2exp_0.4.0          
##  [63] Rcpp_1.0.3                plyr_1.8.4               
##  [65] base64enc_0.1-3           progress_1.2.2           
##  [67] classInt_0.4-2            prettyunits_1.0.2        
##  [69] rpart_4.1-13              deldir_0.1-15            
##  [71] zoo_1.8-6                 ggmap_2.6.2              
##  [73] haven_2.1.1               cluster_2.0.7-1          
##  [75] openxlsx_4.1.0            gmodels_2.18.1           
##  [77] lmtest_0.9-36             hms_0.4.2                
##  [79] mime_0.7                  evaluate_0.14            
##  [81] xtable_1.8-4              XML_3.99-0.3             
##  [83] rio_0.5.10                jpeg_0.1-8               
##  [85] readxl_1.1.0              compiler_3.5.0           
##  [87] maps_3.3.0                KernSmooth_2.23-15       
##  [89] crayon_1.3.4              htmltools_0.3.6          
##  [91] mgcv_1.8-23               later_0.8.0              
##  [93] spdep_1.1-3               Formula_1.2-3            
##  [95] visdat_0.5.3              expm_0.999-4             
##  [97] lubridate_1.7.4           DBI_1.0.0                
##  [99] magic_1.5-9               proj4_1.0-8.1            
## [101] MASS_7.3-51.4             sf_0.8-1                 
## [103] boot_1.3-20               compositions_2.0-0       
## [105] Matrix_1.2-14             car_3.0-0                
## [107] cli_2.0.2                 adegraphics_1.0-15       
## [109] gdata_2.18.0              parallel_3.5.0           
## [111] igraph_1.2.5              pkgconfig_2.0.3          
## [113] rncl_0.8.3                geosphere_1.5-7          
## [115] foreign_0.8-72            laeken_0.5.0             
## [117] sp_1.3-2                  xml2_1.2.2               
## [119] webshot_0.5.1             rvest_0.3.4              
## [121] digest_0.6.25             rmarkdown_1.15           
## [123] cellranger_1.1.0          htmlTable_1.13.1         
## [125] curl_4.1                  shiny_1.3.2              
## [127] gtools_3.8.1              rjson_0.2.20             
## [129] lifecycle_0.2.0           jsonlite_1.6             
## [131] carData_3.0-1             mapproj_1.2.6            
## [133] seqinr_3.6-1              fansi_0.4.1              
## [135] pillar_1.4.2              survival_3.1-11          
## [137] httr_1.4.1                DEoptimR_1.0-8           
## [139] glue_1.3.1                zip_1.0.0                
## [141] png_0.1-7                 class_7.3-14             
## [143] stringi_1.4.6             latticeExtra_0.6-28      
## [145] e1071_1.6-8

Beauchard, O., Verı'ssimo, H., Queirós, A. & Herman, P. (2017). The use of multiple biological traits in marine community ecology and its potential in ecological indicator development. Ecological Indicators, 76, 81–96.

Taugourdeau, S., Villerd, J., Plantureux, S., Huguenin-Elie, O. & Amiaud, B. (2014). Filling the gap in functional trait databases: Use of ecological hypotheses to replace missing data. Ecology and evolution, 4, 944–958.

A. Boyé

12 novembre, 2020